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ABSTRACT. 

A time dependent numerical solution of the linearized continuity 
and momentum equation is developed for sound propagation In a two- 
dimensional straight hard or soft wall duct with a sheared mean flow. 
The time dependent governing acoustic-difference equations and boundary 

i 

conditions are developed along with a numerical determination of the 
maximum stable time increments. The analysis begins with a harmonic 
noise source radiating into a quiescent duct. This explicit iteration 
method then calculates stepwise ?n real time to obtain the transient 
as well as the "steady" state solution of the acoustic field. Example 
calculations are presented for sound propagation in hard and soft wall 
ducts, with no flow and with plug flow. Although the problem with 
sheared flow has been formulated and programmed, sample calculations - 
have not yet been examined. So far, the time dependent finite dif- 
ference analysis has been found to be superior to the steady state 
finite difference and finite element techniques because of shorter 
solution times and the elimination of large matrix storage require- 
ments. 
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INTRODUCTION 


Both finite difference and finite element numerical techniques 
(refs. I to 27) have been developed to study sound propagation with 
axial variations in Mach number, wall impedance, and duct geometry as 
might be encountered in a typical turbojet engine. Generally, the 
numerical solutions have been limited to low frequency sound and short 
ducts, because many grid points or elements were required to resolve 
the axial wavelength of the sound. As shown in reference 1 (eq. (77)) 
for plane wave propagation, the number of axial grid points or ele- 
ments is directly proportional to the sound frequency and duct 1'ength, 
and inversely proportional to the difference of unity minus the Mach 
number (ref. 2). This later dependence also severely limits the ap- 
plication of numerical techniques for high Mach number inlets. 

Customarily, the pressure and acoustic velocities are assumed to 
be simple harmonic functions of time; thus,' the governing linearized 
gas-dynam?c equations (ref. 28-, pg. 5) become independent of time. 

The matrices associated with the. numerical solution to the time inde- 
pendent equations must be solved exactly using such methods as Gauss 
elimination. As a result, large arrays of matrix elements must be 
stored which tax the storage capacity of even the largest computer. 

In an unpublished work at NASA Lewis Research Center by. the author 
using reference 29, as well as in the work of Quinn (ref. 23, pg. 3), 
the matrix has been modified to allow iteration techniques; unfortu- 
nately, the convergence is too slow to be of any practical value. 
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Other approaches, such as in reference 30, might still offer iterative 
possibi 1 ities. 

Some special techniques have been developed to overcome the dif- 
ficulties of low frequency, short ducts, and low Mach numbers.. As 
shown in references 3 and 10, the wave envelope numerical technique 
can reduce the required number of grid points by an order of magnitude. 

In reference 20, this technique was used to optimize multi-element 
liners of long lengths at high frequencies. At the present time, this 
technique has been applied only to the simple cases of no flow and plug 
flow. A numerical spatial marching technique was also developed in 
references 15 and 18. Compared to the standard finite difference or 
finite element boundary value approaches, the numerical marching tech- 
nique is orders of magnitude shorter in computational time and required 
computer storage. The marching technique is limited to high frequencies 
and to cases where reflections are small. 

As an alternative to the previously developed steady state theories, 
a time dependent numerical solution is developed herein for noise propa- 
gation in a two-dimensional soft wall duct with parallel sheared mean 
flow. Advantageously, matrix storage requirements are significantly 
reduced in the time dependent analysis. The analysis begins with a 
noise source radiating Into an initially quiescent duct. This explicit 
method calculates stepwise in real time to obtain the transient as well 
as the "steady" state solution of the acoustic field. The total time 
required for the analysis to calculate the "steady" state acoustic field 
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will determine the usefulness of the time dependent technique. 

Time dependent numerical techniques have been applied extensively 
to one-dimensional sound propagation (ref. 31, pg. 258 and ref. 32), 
two-dimensional vibration problems (ref. 33, P9-; ^52) and the more 
general problem of compressible fluid flow (ref. 3*0* References 35 
and 36 discuss in detail the stability of time dependent numerical solu 
•t'ions. 


In a companion paper to this work (ref. 37), the explicit time 
iteration techniques of references 31 to 36 are extended to include 
sof t-waj 1 impedance boundary conditions which would be encountered 
in inlets and exhaust ducts of turbo-fan jet engines. The analysis 
applies to two-dimensional straight hard and soft wall ducts without 
flow. In the absence of mean flow, the governing acoustic equation 
was the classic second order wave- equation. The time dependent solu- 
tion to the wave equation was found to be superior to the conventional 
steady numerical analysis because of much shorter solution times and 
the elimination of large matrix storage requirements. 

When parallel shear flow occurs in a duct, the second order formu- 
lation of reference 37 cannot.be used’ since the governing differential, 
wave equation is third order (ref. 28, pg. 9). Rather than attempt to 
solve a third order wave equation, in the present paper the first 
order equations of continuity and momentum will be solved simultane- 
ously. First, the governing acoustic equations and boundary conditions 
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are presented for time dependent propagation in a parallel sheared mean 
flow. Next, the governing acoustic-difference equations and boundary 
conditions are derived. The two-dimensional stability theory is then 
used as a guide to estimate the maximum stable time increment. Immedi- 
ately following the mathematical development, numerical solutions are 
presented for one-and-two-dimensional hard and soft-wall ducts. The 
results are compared with the corresponding steady state analytical 
results. Finally, the times required to perform both the time dependent 
and steady state analyses are compared for increasing, number of grid 
points. 
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I. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 


The propagation of sound in' a two-dimensional rectangular duct,- 
shown in Figure 1, Is described by the linearized continuity and mo- 
mentum equations and the appropriate impedance boundary conditions. 


A. Continuity and Momentum 


The linearized equations for mass and momentum conservation (ref. 

I 

28, pg. 5) can be written for a Cartes ianrcoordinate system in the 
following dimensionless form: 
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These and other symbols are defined in the Nomenclature. The dimension- 
less frequency q is defined as: 


_ H « 

” “ 2? ^ 


* * 

H f 
- — — 
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The asterisks define dimensional quantities. 

The foregoing dimensionless equations apply to the scaled Cartesian 
coordinate system in which the dimensionless height (y /H ) ranges be- 

JL JC 

tween 0 and 1 and the dimensionless length ranges between 0 and L /H. . 
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B. Wall Boundary Conditions 

The boundary condition on the transverse acoustic velocity v at 
the surface of a sound absorbent soft-wall duct can* be expressed in terms 
of a specific acoustic impedance Z> 

v(x, 1 ,t) ~ (5) 

where the complex specific wall impedance is defined as 

5 “ V* - 0 + 'X (6) 

p ; c o 

At the lower wall, the sign on Z, is changed to account for the vector 
nature of v. 

In addition to equation (5), another form of the wall impedance con- 
dition will be considered. Substituting equation (5) into equation (3), 
noting that the Mach number is zero at the wall,- and assuming £ is 
a constant independent of time, then 

If. . -alf. (?) 

8y T, 

This form of the impedance boundary condition cannot be used for plug 
flow when a soft wall is present. 

C. Entrance Condition 

The boundary condition at the source plane P(o,y,t) can be of any 
general form with both transverse variations in pressure and multiple 
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frequency content. However, the numerical technique will be compared later 
to previous solutions in which the pressure and acoustic velocities were as- 
sumed to be plane waves at the entrance and to vary as e or in dimen- 

sionless form as e' 27Tt . Therefore, the source boundary condition used here 
is 

P(°,y,t) = e' 27rt (8) 

D. Exit Impedance 


in a manner similar to the wall impedance, the axial acoustic velocity 
at the duct exit can be expressed in terms of a specific acoustic exit- im- 
pedance ? as 


u<L*/H*,y,t) - p < L y.r> t > 


(9) 


For the plane wave propagation to be considered herein, Z, is taken to be 

6 

equal to unity, which- is exact for plane wave propagation in an infinite 

hard-wall duct. Also, choosing £ to be equal to unity has lead to close 

6 

agreement between numerical and analytical results for plane wave propaga- 
tion into a soft-wall duct (refs. I and 3)* More general values for the 
exit impedance can be found in references 7, 15 (eq. B-A) , 16, or 18 (fig. 
7). 


E. Initial Conditions 

For times equal to or less than zero, the duct is assumed quiescent; 
that is, the acoustic pressure and velocities are taken to be zero. For 
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times greater than zero, the application of the noise- source, equation 
(8), will drive the pressures' in the duct. 

F. Complex Notation 

Because of the introduction of complex notation for the noise 
source and wall impedance, all the dependent variables are complex. 

The superscript (i) will represent the real term while (2) will repre- 
sent the imaginary term 

P « + iP^ (10) 

A similar notation applies to the acoustic velocities. 
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II. DIFFERENCE EQUATIONS 


Instead of a continuous solution for pressure in space and time, 
the finite difference approximations will determine the pressure at 
isolated grid points in space as shown in Figure 1 and at discrete time 
steps At. Starting from the known initial conditions at t=0 and the 
boundary conditions, the finite difference algorithm will march-out the 
solution to later times. 

A. Central Region 

Away from the duct boundaries of Figure 1, the derivatives In the 
governing equations can be represented by the following differences in 
time and space: 



k 

v. . 3M 
i >J 

n 3y 
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( 13 ) 



where i and j denote the space indices, k the time index and Ax, Ay, At 
are the space and time mesh spacing respectively. All spacings are as- 
sumed constant. The time is defined as 


k+1 k 

t = t + At o (k + l)At (14) 


Solving equations (11), (12), and (13) for the acoustic pressures and 
velocities yields 


continuity 


k+1 

p . . 
• >J 



J_ At / k 
2q Ax ^ U ?+l,j 



J_ At / k 
2n Ay \ V i,j+1 



(15) 
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2q Ax i+1 , J 



x- moment urn 
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2n Ax 
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( 16 ) 
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y-momentum 


k+1 k / k+1 k+1 \ 

V ',j = V i,j “ 2n^y( P i , j-M “ P i , j~l ) 




(17) 


Equations (15), (16), and (17) are algorithms which permit marching 
out solutions from known values of pressure and velocities at times assoc- 
iated with k and k-1 . First, equation (15) is solved for at all grid 

points prior to solving equations (16) and (17). In this manner, new values 
k+1 

of P are available for use in equations (16) and (17). The procedure is 

explicit since all the past values of P^ are known as the new values of 
k+1 

•P are computed. For the special case at t=0, the values of the pres- 

sure and velocities associated with the k-1 value are zero from the assumed 
initial conditions. 


B. Wall Condition 


Recall that at the wall, the transverse velocity is governed by the 
wall impedance as given by equation (5). Therefore, at the wall, the 
y-momentum equation (17) was initially replaced by the difference form of 
equation (5) which in this case is 


k+1 

v. . 

i ,J 




k+1 

P i,j/C 


( 18 ) 


Thus, equations (15), (16), and (18) are the governing equations at the 
wal 1 . 


The numerical solutions remained stable in the example problems in- 
volving hard-wall ducts (C‘ =co ). Unfortunately, in the soft-wall example 



problems considered, the numerical procedure went, unstable. Introducing, 
equation (18) appears to lead to an instability for finite values of .£. 
An explanation of why this went unstable wi.11 be presented later in the 
section which discusses the stability of the difference equations. In 
the meantime, the following new wall boundary condition is developed to 
replace equation (.18). 

Associated with the new* wall boundary condition,, the grid structure 
shown in Figure 2 is used to establish the pressure gradient at the wall 
boundary which satisfies the impedance condition.. Equation (7) can be 
written in difference form as 

k+1 k 

P. . - P. . 

!_>J . k, i 

At 



Notice that the values of pressure gradient are expressed in the new 

values of pressure. This- can be extremely important in stability consid- 

k+1 

erations, as will be discussed later.- The value of P. . are first found 

» »J 

k+1 

in the central points away from the boundary; consequently, P. is 

a known value in equation (19). Solving for the acoustic pressure at 
the wall yields (j = J) 


k+1 

P '-J 


pk+1 , 1 Ajr p k 

i ,j-l C At t , 


( 20 ) 


A similar equation applies at the lower wall. 
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Equations (20) and (18) now determine the pressure and velocity in 
the wall. The pressures and transverse velocities at the wall would be 
the average of their values at the grid points which straddle the wall. 
The calculation of the axial velocities at the wall grid point is not 
yet required. The axial velocity at the wall will be the value associ- 
ated with the point J-l. 

C. Entrance Condition 

At the entrance conditions, the continuity equation is replaced by 
equation (8) and the difference form' of the momentum equations must be 
expressed in terms of forward differences. In this case, the governing 
difference equations become 



D. Exit Condition 


Similarly to the entrance condition, the governing equations must be 
expressed in terms of the backward differences. In this case, the governing 
difference equations become (i = I) 


k+1 

P. . 
' >J 
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r]Ax 
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(25) 


k+1 


k+1 

v. . 

i ,J 


k+1 

P. ./? 


e 


k At 4 k+l 
v i,j ' 2nAy p,J+t 




E. Spatial Mesh Size 


(26). 


The mesh’ spacing Ax and Ay must be restricted to smal.l values to re- 
duce the truncation error.- To resolve the oscillatory nature of the pres- 
sure, the required number of grid points in the axial direction suggested 
is (see ref. 1) 


12qL 
~(1 + M) 


(27) 


while the number of points in the transverse direction suggested in reference 
37 is 

J > 12n (28) 
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III. STABILITY 


In the explicit time marching approach used here, round-off errors 

<• 

can grow in an unbounded fashion and destroy the solution If the differ- 
ence equation is improperly formulated or the time increment At is taken 
too large. Although at least six methods (ref. 3*1, pg. 48) of stability 
analysis exist, none are entirely adequate when obtaining actual solutions 
to differential equations (ref. 3*t> pg. 51) • Consequently, numerical ex- 
perimentation will be used to determine the actual stability. The simpler 
stability analyses will be used to guide the development of the differ- 
ence form of the governing equations and the choice of the time increments. 

* 

Using the two-dimensional (space and time) no flow stability analysis 
of Courant, Friedrichs and Lewy (ref. 31 > pg- 262) as a guide, p^*^ values 
were used in equations (12) and (13) instead of P^. If the values of P 
at pk instead of at P^ had been used in equations (12) and (13), the 
iteration scheme would be unstable. To avoid instabilities, as a general 
rule, the new values of the dependent variables will be used whenever pos- 
sible' (see equation (19), for example). 



At < 
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(30) 



To account for Mach number, equation (30) is empirically modified to 



where the largest value of the local Mach number is used. Equation (31) 
is used to set the initial time increment. If an instability should occur.,- 
the time increment will be reduced until stability occurs. 


20 



IV. STEADY STATE PRESSURES 


In the sample problems to be presented in the next section, the time 
dependent results will be compared to the results of the steady harmonic 
solutions of reference 10, The purpose of this section is to show the 
rationale for constructing a steady state solution from the time dependent 
resul ts. 


A. Steady Harmonic Solution 


The steady harmonic pressure p g (x,y) is-defined as a solution to 
equations (1) to (3) when the pressure' is assumed to' be a simple harmonic- 
function of time: 


P (x,y,’t) = p s (x,y)e ,2lTfc (32) 

In this case where the source is a- s-imple harmonic function of time, p g 
represents the Fourier transform of p(x,y,t) (ref. 28, pg. 11). 


For a semi-infinite duct (or an equivalent finite duct with a- p"c" 
exit impedance) with plane wave propagation, uniform Mach number, and hard 
walls,, the solution for p g is (ref-. 10, eq. (18)) 

“ i 2717] X 


p 


s 


e lHHT 


(33) ■ 


In the next section, a transient solution to this problem will be compared 
to equation (33) . 
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B. 


Transient Times 


Recall, at the start of the numerical calculations, the acoustic 

pressures and velocities were assumed zero throughout the duct. Then, 

a pressure source begins a harmonic oscillation at x^O for t greater 

than zero. For the special case of a plane wave propagation without 

flow in a hard wall semi-infinite duct, the analytical solution to the 

wave equation indicates that the transient solution ends and the steady 

solution p begins when (ref. 38 , pg. 305) 
s 

t t > T|X (34) 

or in terms of real variable when 

JL JL JL 

t? > x7c~ (35) 

t o 

The transient time in equations (3*0 and (35) represents the time for 

JL JL 

the wave to travel down to the end of the duct, x = l\ 


The transient time can be shortened or expanded depending on the 
direction- of the flow and the transverse velocity distribution, in 
general, the length of time for the transient to occur will be found by 
a trial numerical procedure. For the simple case of plug flow, equation 
(3*0 is now modified to include a uniform Mach number. 


t 


t 


> 


nx 

i + M 


( 36 ) 


Therefore, for the special' case of one ^dimensional plane wave propagation, 
the initial transient is assumed to' pass when equation ( 36 ) holds. As a 
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factor of- safety In the present calculations, the transient calculation 
will be continued into, the steady domain for one period of oscillation 
before the Fourier pressure p is calculated.. Therefore, in this paper 


t = 
c 


JL JL 


til'Vh" 

Tl + M) 


+ 1 


(37) 


and 


p(x,y). « 


P(x,y,t) 

I2irt c 
e c 


( 38 ) 


For more complicated problems, such as with higher order modes, 
where reflections are important, or where complicated flow gradients 
exist, t should be increased in successive steps to check for conver- 

Q 

gence. 
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V. SAMPLE CALCULATIONS 


In the sample problems to follow, the time dependent results will 
be compared to the analytical results of equation (33) and the steady 
harmonic numerical solutions of reference (10). 

A. Hard-Wall Duct 

Numerical and analytical values of the pressure p(x,y) are computed 
for the case of a hard-wall duct for plane wave propagation with £ « 1 

exit impedance (equivalent to a semi-infinite duct) for no flow and plug 
flow with Mach numbers of -0.5 and +0.5. The calculations are made with 
a length to height ratio (L“/H) of 1 and a dimensionless frequency ri 
of 1. The analytical and numerical values of the acoustic pressure pro-; 
files along the duct are shown in Figures 3, A , and 5. As seen in these 
figures, the agreement between the analytical and numerical' theory is 
reasonably good. 

Some inaccuracy exists in the pressure at the entrance (x = 0). As 
seen in Figures 3 and 4, the pressure at the second grid location takes 
a slight jump. This is believed to result because the approximations 
for the spatial derivatives are only first order at the corners. Future 
work should be concerned with developing higher order difference approxi- 
mations which are more accurate and stable. Higher order difference 
would also be desirable in order to increase the mesh spacing used in 
these calculations, and thereby reduce the number of grid points. 



B. Soft-Wall Ducts 


As another example of the time dependent analysis, the pressure 
distribution was computed for the case of plane wave propagation with- 
out steady flow and with a 5^ » 1 exit impedance and a wal 1 with 
impedance value of 0.16 - i 0.3^. The calculation was made with a 
length to height ratio of 0,5 and a dimensionless frequency of 0.6. 

The results of the time dependent analysis along with the results of 
the solution of the equivalent steady state Helmholtz equation are dis- 
played in Figure 6. The numerical results for the steady spatial solu- 
tion P s (x,y) are tabulated in Appendix F of reference 10. 

Again, as seen in Figure 6, the steady state and time dependent 
solutions are in reasonable agreement. Hopefully, the difference be- 
tween both theories- can be resolved by using higher order difference 
approximations. 

C. Grid Point Variations 

Figure 7 shows the effect of increasing the number of grid points 
on the computational time of the time dependent approach for the hard 
wall duct associated with Figure 3. Roughly, as seen in Figure 4, the 
computational time is proportional to the number of grid points used. 
This is considerable advantage over the steady state technique in which 
the computational time more nearly increases with the square of the 
total grid points. 
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Figure 8 shows a comparison of the computational times for the 
steady state and time dependent problems associated with the hard' wall, 
duct shown in Figure 2. For this no-flow example, the continuity and 
momentum equations, (1) to (3), can be combined, into a single wave 
equation 


2 <£l 

9t 2 


T1 




(39) 


The time dependent solution to equation (39) was presented in reference, 
37 and is represented by the lowest line in Figure 8„- 

As seen In Figure 8 for J = 20, the time dependent analysis pre- 
sented herein is roughly equal to the steady state analysis. The value 
of J (grid points in y-di rection) was restricted to 20 because of prac- 
tical limitations on the size of the matrix which could be effectively 
handled in the steady state analysis. When J was increased to 50 wi.th 
50 axial grid points, the steady state analysis required 5500 seconds as 
compared to less than 150 seconds for the time dependent analysis. This 
large increase- in the steady state solution time results because of the 
manner in which the general matrix was partitioned (ref. 10, pg. 1A). 

In this case, the storage and computational times are proportional to the 
total number of transverse grid points squared. 


As seen tn Figure 8, the time dependent solutions for the continuity 
and momentum formulation in this paper require computational times a fac- 
tor of 10 greater than the transient solution times to the wave equation 
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for the same degree of accuracy. This results- because much smaller time 
increments At must be taken with the continuity and momentum solutions 
to obtain the same accuracy or truncation error. In the numerical formu- 
lation of the wave equation {eq. (39)), the second derivative of time is 
expressed in terms of the usual central difference approximation: 



p!? + ! 

i .J 


2P^ . + P? ! 
Li 1 

-.2 




0(At 2 ) 


(40)- 


2 

where the truncation error is of order At . On the other hand, in the 
continuity equation (eq. (11)), the first derivative of time is expressed 
as a forward difference 


9P 

9t 


-- pk +1 _ p k 
_ _ r • » r » • 

= ‘»J '»J 

At 


+ 0(At) 


m 


where the truncation error is of order At. Consequently, a numerical 
solution based on equation (Al) will require smaller time steps for the 
same degree of accuracy. 

D. Shear Flow 


Although the shear flow difference equations have been programmed, 
at the present time shear flow example problems have not been examined. 
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CONCLUSIONS 


A time dependent two-dimensional explicit numerical procedure was- 
developed for the parallel sheared mean flow form of the separate con- 
tinuity and momentum equations. This time marching technique was .found 
to be stable for both no flow and plug flow. At the present time, shear 
flow examples have not been attempted. A special wall boundary condi- 
tion was developed to insure stability for the soft wall case. 

With the possible exception of the wave envelope technique (ref. 

10) or the spatial- marching technique (ref. .18), the numerical time de- 
pendent method of analysis represents- a significant advance over previous 
steady numerical theories. By eliminating large matrix storage require- 
ments, numerical calculations of high sound frequencies are now possible. 
Because manipulation of matrices is omitted, the time dependent approach 
is much easier to program and debug. 
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Figure L - Grid-point representation of two-dimensional flow 
duct 



Figure 2. - Grid-point representation of two- 
dimensional flow duct using pressure boundary 
condition at walls. 
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Figure 3. - Analytical and numerical pressure profiles 
for one-dimensional plane wave sound propagation 
in a hard wall duct for M c 0, rj = 1, and HH* = 1 
(I - 25, J - 20, At = a 004084). 
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Figure 4 - Analytical and numerical pressure profiles 
for one-dimensional plane wave sound propagation 
in a hard wall duct for M = +05, r? = 1, and 
L' : 7Ff = 1 (I = 25, J = 10, At - Q 002438). 
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Figure 5. - Analytical and numerical pressure profiles 
for one-dimensional plane wave sound propagation 
in a hard wall duct for M = -0.5, n = 1, and 
L*/H* = 1 (1 « 49, J » 10, At = a 00064). 
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Figure 6. - Pressure magnitude on the duct axis for incident 
plane wave sound propagation in a soft wall duct (M = 0, 
r] - 0.6, V'lH" - .5, Z = 0. 16 - id 34, 1 = 6, J - 10, 

At ■ . 0056). 
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Figure 7. - Effect of increasing number of grid points 
on calculational time of transient solution for plane 
wave propagation in a hard wall duct (M E 0, ■/} = 1, 
L*/H* - 1, and a = 1). 
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Figure 8. - Comparison of transient and steady state 
calculations for plane wave propagation in a hard wall 
duct <M = 0, v = 1. L*/H* - 1, J - 20). 
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